Corrections to Scaling in the Hydrodynamic Properties of Dilute Polymer Solutions 



Burkhard Diinweg*, Dirk Reith, Martin Steinhauser, and Kurt Kremer 
Max-Planck-Institut fiir Polymerforschung 
Ackermannweg 10, D-55128 Mainz, Germany 
(February 1, 2008) 

We discuss the hydrodynamic radius Rh of polymer chains in good solvent, and show that the 
leading order correction to the asymptotic law Rh oc N" (N degree of polymerization, v « 0.59) 
is an "analytic" term of order A^ - ' 1- "- 1 , which is directly related to the discretization of the chain 
into a finite number of beads. This result is further corroborated by exact calculations for Gaussian 
chains, and extensive numerical simulations of different models of good-solvent chains, where we find 
a value of 1.591 ± 0.007 for the asymptotic universal ratio Rc/ Rh, Rg being the chain's gyration 
radius. For chains the data apparently extrapolate to Rg/Rh ~ 1.44, which is different from 
the Gaussian value 1.5045, but in accordance with previous simulations. We also show that the 
experimentally observed deviations of the initial decay rate in dynamic light scattering from the 
asymptotic Benmouna-Akcasu value can partly be understood by similar arguments. 



I. INTRODUCTION AND SUMMARY 

It is well-known that the average size R of an iso- 
lated flexible uncharged polymer chain in good solvent is 
asymptotically proportional to N v , where N is the de- 
gree of polymerization, and v » 0.5877EJ. This law holds 
for any measure of the chain size, the most popular of 
which are the mean square end-to-end distance, 



(Rl) = (rl N ) , 
the mean square radius of gyration, 
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ij 

is measured. In the small q limit, it decays like 
S{q,t) 
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where D is the chain diffusion constant, .which, within 
the framework of Kirkwood-Zimm theoryo is related to 
Rh via 



N 6707 \ Rh / 
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and the hydrodynamic radius 



.v-E\,. 



Rh I 
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In these equations, we have assumed that the chain is 
composed of N monomers (i. e. N— 1 bonds) at positions 
r^f, 4 = 1,..., JV, and rij = |rj — fj\. Experimentally, the 
gyration radius is determined from small-angle scattering 
experiments; for small wave numbers q the single-chain 
static structure factor behaves likeB 



S ( q } = ^E^ exp W'^ij]) 



N 
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Conversely, the hydrodynamic radius is determined via 
small-angle dynamic light scattering experiments, where 
the dynamic structure factor 



where Dq is the monomer diffusion constant (usually this 
contribution is neglected), ks is Boltzmann's constant, T 
the temperature and 77 the solvent viscosity. In principle, 
it is possible to obtain Rh, as a static quantity, also from 
purely static scattering, as is seen from the relation 



\Rh/ 



dq(S(q)-l) 



(8) 



but for technical reasons, this has so far not been applied 
in experiments. 

When analyzing data for the chain size, one has to 
take into account that the law R oc N v holds only in the 
asymptotic limit N — > 00, while for finite chain lengths 
deviations occur. This is particularly important for com- 
puter simulations, where data with high statistical accu- 
racy can be obtained. For this reason, corrections to scal- 
ing have been worked out in great detail, and exploited 
in high-resolution numerical studies, for the end-to-end 
distance and the gyration radius, where the relation 



"corresponding author, Electronic Mail: duenweg@mpip-mainz.mpg.de 



1 



(R%) = AN 2 " 1 



B 



(9) 



(and analogously for Rq) holdsa. Here A and B are non- 
universal amplitudes, while A is a universal correction- 
to-scaling exponent, whose value is difficult to determine 
beyond the accuracy A ps 0.5 (A pa 0.56 according to Ref. 
|l|, A w 0.43 according to Ref. |]). The omitted terms are 
further powers N~ A ^, A^ A_2 , . . . , as well as N~ A2 , 
2_1 , . . . (i. e. there are further larger correction-to- 
scaling ejsponents) , plus so-called "analytic" terms -/V -1 , 
7V~ 2 , . . .a. The important point to notice is that the "an- 
alytic" corrections will arise even for a Gaussian chain, 
and are due to the fact that the chain consists of a finite 
number of b eads. This will be demonstrated explicitly 
in Sec. II A. Conversely, the "non-analytic" corrections 



are due to the fact that, in the language of renormaliza- 
tion group (RG) theory, the chain's Hamiltonian is not 
identical to the fixed point Hamiltonian. The exponent 
A is related to the largest sub-leading eigenvalue of the 
RG transformationj-at the fixed point. In first order e ex- 
pansion its value istl A = vu> with \jv = 2 — e/4 + 0(e 2 ) 
and to = e + 0(e 2 ), where e = 4 — d aad d is the spatial 
dimension. Higher-order calculationsu have resulted in 
A = vuu = 0.588 x 0.82 = 0.482. We adopt here the con- 
vention (which we view as quite natural) to distinguish 
the terms by their different origins, and call correction 
terms "analytic" corrections if they are present even in 
the Gaussian limit, while we call "non-analytic" correc- 
tions those terms which occur exclusively for excluded- 
volume chains. As we will see in Sec. II B, "analytic" 
corrections defined in this way do not necessarily imply 
integer powers of N . 

As for the hydrodynamic radius of self-avoiding walks 
(SAWs), there is no high-resolution numerical study 
available, and corrections to scaling have not yet been 
dealt with systematically. This is somewhat unfortunate, 
since the corrections are unusually large for Rh, and of 
experimental relevance. For a good solvent chain, one ex- 
pects again N~~ A etc. terms, plus "analytic" corrections. 
It is the main purpose of the present paper to show that 
the leading-order term of these latter corrections is now 
given by 
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whe re B is usually positive. We will show in Sees. II B 
and II C that this form is a straightforward consequence 
of discretizing the chain into beads. As 1 — v « 0.41, 
this will, for long chains, ultimately dominate over the 
A^~ A term, where the exponent is (according to Refs. [j], 
^ and ^) slightly larger. Nevertheless, the exponents are 
so close that in most cases one will observe contributions 
from both terms. On the other hand, it is a well-known 
empirical fact that for many experimental systems, as 
well as for most computer models, the corrections to scal- 
ing of the gyration radius are quite weak, such that the 



A^~ A term should have a rather small amplitude. One 
could therefore expect that the corresponding amplitude 
of the A" _A contribution in Rh is also quite small. Then 
the most likely candidate for explaining the experimental 
and numerical observation that Rh is usually subject to 
very large corrections to scalingcl would actually be the 
"analytic" AM 1- ") term. 
For a Gaussian chain we are able to solve the problem 



exactly, see Sec. II A 
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N~ 1/2 ( l-BN~ 1/2 
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with B = -(3/4)C(l/2) » 1.095266 (here ( is Riemann's 
zeta function), and b denoting the root mean square bond 
length. As v = 1/2 for a Gaussian chain, this form is con- 
sistent with Eq. [l(]. 

The difficulties in observing the asymptotic N ''.-scaling 
of Rh have a long history. Adam and DelsantiuB per- 
formed dynamic light scattering experiments and found 
an effective power law Rh oc TV 0,55 . This is quite typical, 
and has been found in many-other experiments, tooolill, 
although an exponent of 0.6lli3 has been reported as well. 
A reduction of the effective exponent is indeed expected, 
as seen from Eq. [l^, rand is also observed in Brownian 
Dynamics simulations!!^. As a caveat, note that a scat- 
tering experiment does not measure Rh, but rather the 
diffusion constant D. This quantity has an additional 
Dq/N contribution (Eq. ^), which is of the same order 
as the leading correction of Eq. |l^. Therefore, the cor- 
rections in D are weaker than those in Rh- Nevertheless, 
the Dq/N term is typically not large enough to fully com- 
pensate the corrections in Rh ■ This is easily seen for the 
Gaussian case from Eqs. [?] and [H]: The monomer diffu- 
sion constant Dq can be written as Dq — kBT/(6wqa), 
which defines a monomer Stokes radius a. Thus 



D 



N- 1 

| (3.6853N- 1/2 - 4.0364AT 1 
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Since a should be of the order of the bond length b, one 
sees that a large A?^ 1 contribution remains. 

A first attempt to explain the experimental observation 
is due to Weill and des Cloizeauxtj. They conjectured 
that u e ff = 0.55 is due to non-perfect solvent quality, 
and a crossover between good solvent behavior (y « 0.6, 
large length scales) and 8 solvent (y = 0.5, small length 
scales). In particular, they pointed out that the averag- 
ing over 1/r assigns a very large statistical weight to the 
small distances. Although this latter argument is true, 
and generally accepted as the basic origin for the slow 
convergence of Rh, the explanation in terms of solvent 
quality has turned out to be incorrect. In Ref. [l5| it was 
demonstrated that Rh should not be much more sus- 
ceptible to solvent quality effects than Re or Rq — the 
enhanced sensitivity of Rh to the small distances is bal- 
anced by the fact that Re and Rq are more sensitive to 
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the decreased swelling of the chain near its ends: A SAW 
is inhomogeneous, i. e. (fy) depends on the position of 
the ij bond on the chain, and is systematically i^fgpr in 
the interior -as-ias been shown both numericallyMEj and 
analyticallyB'ta _ 

Furthermore, Schafer and BaumgartnerEJ performed a 
detailed RG calculation and predicted in one-loop order 
for the universal amplitude ratio 



lim 



Rg(N) 



1.06 x 



30? 



1.06 x 1.5045 « 1.595; 



(13) 



here 8/(3y / 7r) is the exact random walkp-CRW) value. 
Other RG studies resulted in = 1.56213 (this value 
was later revised to 1.51, see Ref. |l|) and p^ — 1.62E3, 
while a semi-empirical relation based on fitting the dis- 
tribution function of internal distances to light scatter- 
ing data yields poo — 1.5955li3ta. A value of p^ ~ L-6 
was also found in Brownian Dynamics simulationslij. 
From Eq. |l^ it is clear that p should be subject to an 
jy-Ci-") correction for finite chain length; nevertheless, 
experiments have so far not reported a systematic de- 
pendence.^ molecular weight. Typically, values around. 
p » 1.5BE3, or p « 1.6 / p « 1.3 for different solventsE3 
are found in the good-solvent regime. In view of the in- 
accuracies of the experiments (Rq typically has an error 
of 59cEj) the inability to observe a systematic behavior 
in N is not very surprising. 

In order to contribute to the resolution of these ques- 
tions, we have performed computer simulations of very 
different models of polymer chains, both for SAWs and 
for O chains, and calculated Rq and Rh, as outlined in 



Sec. Ill . To provide a complete and well-converged data 
set represents the second main goal of our paper. To 
our knowledge, our results are the most accurate data 
obtained on Rh so far. Concerning Rq, however, our 
data are less accurate than those of Li et alu, which 
we still view as the most precise numerical study on the 
SAW problem so far. Therefore we have taken their val- 
ues for the exponents v and A for our fits. We find 
Poo = 1.591 ±0.007 for good-solvent chains, in very good 
agreement with Ref. [TBI. (Note that our error estimate is 
probably overly optimistic, since it only includes statis- 
tical errors and completely neglects systematic errors.) 

Theoretical and numerical investigations on correc- 
tions to scaling in Rh have first iacpsed on the RW 
case. The work by Guttman et alSSc3 showed by an- 
alytical calculation that a Gaussian chain should— obey 
Eq. [fl|. The prefactor of the correction was firstE -de- 
termined only approximately, B « 1.125, while latere 2 ] it 
was given exactly in terms of an integral. Furthermore, 
Monte Carlo (MC) simulations of lattice chains at the 
O point revealed that in this case the ratio Rq/Rh ap- 
parently does not converge to its Gaussian value 1.5045, 
but rather to roughly 1.4. Our simulations (see Sec. Ill) 
find a similar behavior (poo w 1.44). We believe that this 



can be explained qualitatively from RG arguments^ as 
follows: The asymptotic behavior is expected to be gov- 
erned by the Gaussian fixed point, and thus poo, as a uni- 
versal amplitude ratio, is expected to assume the Gaus- 
sian value. However, the numerical extrapolation will 
only produce this value if all relevant correction terms, 
i. e. the analytic A" 1 / 2 term, plus the non-analytic cor- 
rections, are consistently taken intp account. Neither the 
data analysis by Guttman et a/BlJO, nor ours, fulfill this 
requirement, as both just fit to Eq. [h] with v = 1/2, 
and thus are expected to produce substantial systematic 
errors in p^ . To do this in a better way is practically im- 
possible, since (i) our O data have insufficient statistical 
accuracy to allow for additional fit parameters, (ii) the 
precise form of the correction terms is unknown for Rh 
(in contrast to Re and Rq, for which the leading-order 
terms have been calculated by tricritical field theoryEj, 
with the interesting feature that they are universal), and 
(iii) the non-analytic corrections vary extremely slowly 
(logarithmically) with N, such that either one would need 
unrealistically long chains to ensure dominance of the 
leading orders, or an expansion up to unrealistically high 
order. These problems have-, been elucidated in quite 
some detail for Re and Rj^, explaining previous dif- 
ficulties in the. interpretation of highly accurate MC data 
on O chainsEJ. Ip,th|is context, it should be mentioned 
that experimentsr 5 H 2 1 typically find a value of p = 1.3, 
i. e. a similar reduction as in the good solvent case. 

Later, MC data were taken of excluded-volume (EVl 
chains with SAW statistics. Schafer and Baurngartnerli3 
used chains of up to 161 monomers, with an EV strength 
particularly close to the SAW fixed point, such that poor- 
solvent effects can be ruled out. The inhomogeneous 
swelling was demonstrated, and the Rh data were fitted 
with Eq. [To]. This was done with an empirical correction- 
to-scaling exponent of 1 /2 instead of 1 — v. The same ev- 
idence was shown in the simulation data by Batoulis and 
Kremei£3 of chains of length of up to N 400. Ladd and 
FrenkelES simulated chains of length of up to N = 1025 
and were able to describe their Rh data via Eq. [l(J with 
A = 3.84 and B = 1.06, but without detailed justifica- 
tion of their use of-,thc correct 1 — v exponent. Schafer 
and Baumgartnerlij concluded from both their analyti- 
cal studies and their simulation data that not the solvent 
quality, but rather the chain's microstructure is res ponsi - 
ble f or the slow convergence. Our reasoning (Sees. 
II C ) , which is similar to the one by Guttman et al\ 
exactly supports this picture: The corrections are due to 
the fact that the chain is discretized into beads, or, in 
other words, to the fact that there is a lower length scale 
cutoff for the frictional properties. However, the notion of 
"stiffness" , which is often used in this contextll3, i s, in our 
view, somewhat misleading: As outlined in Sec. II D. we 



expect a large local chain stiffness to decrease the correc- 
tion until it ultimately even changes its sign. The .same 
conclusion has been found by Akcasu and Guttmant 2 ] for 
stiff chains without excluded volume. 
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In the context of dynamic light scattering of dilute 
polymer solutions there is yet another unresolved puz- 
zle. As Akcasu et al. have shownBS, the initial decay 
rate of the dynamic structure factor, 



can be written as 



d S(q,t) 



dt S(q, 0) 



t=o 



(14) 



and 



(4> = 6 2 N-i| 



(ry 1 > = 6 1 /V- 1 /»6- 1 |i-jT 1/a : 



and hence 



(18) 



(19) 



Q(q) 



Ejj \<t- D ij ■ qexp(iq- ' 

E» ( ex p(«f-^j)) 



(15) 



(R 2 E )=b 2 (N-l) = b 2 N (l-l) 



(20) 



where is the diffusion tensor. Equation [T^ is a an( j 
rigorous result, the only assumption being that the 
chain dynamics can be described by Kirkwood's diffusion 
equationi. Usually, Dij is taken as the Oseen tensor, 



D^ =D Sijl + (!-<%) 



k B T 



(i+hj®hj), (16) 



where fi 



denotes the tensor product of the unit 



vector in fij direction with itself. In this case, Eq. |15| is 
just thjj-o > generalization of Eq. |t]. It can then be 
shown0G!J that for q in the scaling regime Rq 1 <C q <C b^ 1 
(b denoting the bond length), or, strictly spoken, in the 
limit qb — ► 0, qRa — * oo, the relation 



(17) 



holds, where the numerical constant C only depends on 
chain statistics: C — 0.0625 for RW statistics (y = 1/2) 
and C = 0.0788 for SAWs (y = 0.6). This has been 
tested hw—Jiftht-. scattering experiment^ both for good 
solventsa'Ed'ElrO and for 9 solventscSO. In both cases 
the relation is verified with reasonable accuracy, but with 
a prefactor C which is systematically smaller than the 
theoretical prediction. The reasons for this shift are not 
clear; an attempt by a generalized theory which intro- 
duces drainingp.3 so far had only limited successful. In Sec. 
HE we show that the deviation can partly be explained 
by the fact that in reality neither qb = nor qRc = oo 
holds. Taking these nonidealities crudely into account, 
we find a shift in the same direction, which is however 
smaller than the experimental one. Nevertheless, we be- 
lieve that this third main result is of direct relevance for 
the analysis of experimental data. There are alsp-.some 
indications from Molecular Dynamics simulations^] that 
the description in terms of the Kirkwood theory is insuf- 
ficient on these length and time scales. 



N-l 



N 2 

= h 2 N (l - -1 
6 V N 2 



71=1 



(21) 



where we have used elementary summation formulae. For 
the hydrodynamic radius, we find analogously 



W> = V^£»- 1/2 ("-»>- 



N-l 



irbN 2 



(22) 



Accordi ng t o the Euler-Maclaurin formula (see Appendix 
0), Eq. [A|, the sums can be expanded as 



Y^n- 1 ' 2 = 2N 1 ' 2 - \n-^ 2 + C(\)+ 0(N~ 3 / 2 ), 



^ 3 2 

n—1 

+ C(-1)+0(N-^ 2 ). 



Hence, 



N-l 



(23) 



Y n" 1/2 (A - n) = \n^ 2 + AC ( \ ) + 0(N°) (24) 



and 



II. ANALYTICAL THEORY 

A. Hydrodynamic Radius of a Gaussian Chain 

For a Gaussian chain with root mean square bond 
length b, we have 



(Rh 1 ) = 



6_8_ 

3b 



N -l/2 



( i+ i c G)^ i/2+o( ^ 3/2) )' (25) 



which is the result anticipated in Eq. |Tl|. 
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B. Hydrodynamic Radius of a Good Solvent Chain 

For a linear SAW, the main dimcultAjrisr-fche fact that, 
unlike for a RW, (rfj) and (r^ 1 ) do noflli§li3 just depend 
on \i — j\, but rather on the positions relative to the ends 
as well. In order to obtain the leading-order analytic cor- 
rections due to discretization, we can restrict the discus- 
sion to the leading-order scale-invariant behavior, i. e. 
we can assume that the SAW is strictly scale invariant 
with the exponent v, with no non-analytic corrections. 
If we would include the latter, they would just generate 
further additive terms in our expressions. In what fol- 
lows, we therefore omit them, for the sake of simplified 
notation, but keep in mind that they have to be added at 
the end in order to obtain the full expressions. We thus 
assume the relations 



100 



i>G (Ax, Xy) 
b H (Ace, Ay) 



\ 2v <j) G (x,y) , 



where we have introduced the notation 



<t>a{i,0) = (rfj) 



(26) 
(27) 



(28) 
(29) 



The definitions of Rg and Rh lead us to study the sum 



N n-1 



(30) 



for <p = <pQ and (f> — (pu, respectively, by means of the 
Euler-Maclaurin expansion of Appendix ^| Treating the 



inner sum first, we find 



4>{m, n) — const. + tp(n) 



(31) 



m— 1 

with the formal expansion 

f n Id 

ip(n) = / dx(j){x,n) + — —(j>(x,n) 



(32) 



since <j>(n, n) vanishes. Note also that the constant in Eq. 
pl| does not depend on n, hence 



a(N) = (N- l)const. + ip(n) 

n=2 
rN+1 

= (N — l)const. + / dyip(y) 



1 



ip(N + 1) + const. + 



(33) 



a; 




FIG. 1. S(q) as defined in Eq. ||, for a chain of length 
N — 64, for model C (see Sec. ill), demonstrating the q~ 1/V 



decay, followed by oscillations around unity. The dashed line 
is the simplified structure factor of Eq. |34]. 



a(N) 



N+l 

dy / dx(j)(x,y) 
2 Ji 

f N+1 1 d 

J 2 dv Y2Tx^ v) 



1 



N+l 



dxcj){x, N + l) 



TVconst. 
const. + 



(34) 



Inserting Eq. p3, we find 



After transformation to the reduced variables u — x/N 
and v — y/N, and exploiting the scaling behavior of </>, it 
is possible to determine the order of each term. For the 
gyration radius, we find 0{N 2+2 »), 0{N 2u ), 0(N 1+2u ), 
0(N 1 ), 0(N°), respectively, for the five terms in or- 
der of their appearance. Conversely, for (R^ 1 ) the dif- 
ferent scaling behavior of cj> implies 0(N 2 ~ lJ ), 0(N~ lJ ), 
OiN 1 -"), OiN 1 ), O{N ) for the corresponding orders. 
For (R G ) the leading order is 0(N 2+2u ), while the next 
sub-leading order is 0(N 1+2l/ ), resulting in a leading or- 
der correction of 0(N^ 1 ). For (i?^ 1 ), the leading order 
is 0{N 2 - U ), followed by the O^N 1 ) term. Thus the cor- 
rection to scaling for Rg has order 0{N~^ 1 ~ V ^). This 
proves Eq. Of course, this consideration does not 

prove that the amplitude B in Eq. [lO] is positive; how- 
ever, this is expected from the result for Gaussian chains. 



C. Alternative Derivation 

Equation [l^ can also be derived in a more heuristic 
way, starting from Eq. ||. Figure |l| shows the typical be- 
havior of S(q): For wave numbers q with R^ 1 <C q <C b^ 1 
the structure factor exhibits a power-law decay q~ x l v 
which indicates the chain's fractal geometry, while for 
larger q it oscillates around unity. We therefore can in- 
troduce a cutoff wavenumber go from which on there is no 
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further contribution to the integral, i. e. go is the small- 
est of all the g's with the property J?° dq(S(q) — 1) = 0. 
Hence, 



1 

Rh~ 



2 

ttAT 



/ dqS{q) - 
Jo 



2go 
ttN' 



(35) 



It is physically clear that for a flexible chain go must be 
roughly (2n)/b, apart from a numerical prefactor of order 
unity. Moreover, the fractal q~ x l u decay of S(q) roughly 
extends up to go, at which point S(q) « 1 is reached. We 
now introduce a modified structure factor S(q), which is 
identical to S(q) up to g = go, but extends the g~ 1//,y 
decay up to g = oo. In this latter regime, we have 



S(q) 



-X/u 



(36) 



where a is a numerical prefactor of order unity. Therefore 
we can write 

4r = ^T d«S(q) -A-rT ddS(q) ~ ^ (37) 



R H nN 



nN 



nN 



Evaluating the second integral, and writing S(q) in scal- 
ing form, 



S(q) = Ns(qR G ) 



(38) 



(here we have again assumed strict scale invariance, i. e. 
absence of non-analytic corrections to scaling, for the 
same reason as outlined at the beginning of the previous 
subsection), one finds 



Rg 2f 2 / v 

= — I axsyx) a 

R H 7t J 7r V 1 - v 



M^, 09) 



i. e. again a negative correction of order 0(N ^ 



D. Effect of Chain Stiffness 



We now assume 



S(q)=a^(± 
9i \qi 



-X/u 



for q > gi, and 



S(q) = p 



(41) 



(42) 



for gi < g < go . Here, a and (3 denote prefactors of order 
unity, and the g _1 decay results from the local stretch- 
ing. Evaluating the integrals, and using Eq. |38|, one thus 
finds 



Rg 
Rh 



dxs(x 
'o 

2 ( 
-— a— 

7T \ 1 



(43) 



v 



v 



1-0^)^. 
qi) N 



In order to compare with the flexible case, we still have 
to take into account that stiffness tends to increase the 
gyration radius, by roughly a factor of (go/gi) 1_I/ : 



Rg 
~Rh~ 



dxs(x) 



(44) 



1-1/ 



+ 1-/3 In 



go 



N ' 



where Rq denotes the gyration radius in the flexible 
case. 

The prefactor of the correction term hence depends on 
the stiffness parameter go/gi in a non-trivial way; for 
small go/gi both an increase and a decrease are possible, 
depending on the parameters. For sufficiently large stiff- 
ness one always obtains a decrease of the correction, and 
ultimately even a change of its sign. 



The advantage of the approach of the previous sub- 
section is that it can be easily generalized to study the 
influence of local structure, since it is well-known how 
this is reflected in S(q). For a locally stiff chain with a 
persistence length large compared to the bond length b, 
one expects that go is roughly unchanged with respect 
to the flexible case. However, the q~ x l v decay does no 
longer extend down to g w go, but only to g w gi, where 
gi is a crossover wave number, whose inverse is a typical 
length scale below which stiffness effects are important. 
With 5(g) being again the continuation of the q~ x l v de- 
cay up to q — oo, we have 



1 

Rh 


ttN , 


/•OO 

/ dqS(q) - 
Jo 


2 

tt/V 




+ ^/V , 


rQo 

/ dqS(q) - 
'1i 


2g 
■kN 



p oo 

dqS(q) 



(40) 



E. Initial Decay Rate 

In this subsection, we are concerned with the initial 
decay rate 0(g), see Eq. |ll| Splitting the sum in the nu- 
merator into diagonal and off-diagonal terms, one finds 



n(g) 



Dog 2 

S(q) 

^£(g-.5,.gexp W . 



(45) 



Following Refs. |,|37|, we use the Fourier representation 
of the Oseen tensor for the off-diagonal elements, 



- _ k B T 1 
ij ~ 77 (2tt) 3 



, 1 — k ® k - 
d k — exp(ifc • nj), (46) 
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to find 



fi(g) 



D q 2 



1 k R T 1 



S(q) S(q) rj {2nf 



d*k 



(k ■ q) 2 



k 2 



(S(k + q) 



- 1 



We now focus attention on the dimensionless quantity 



C(q) = 



n 



q 3 k B T 



Q(q) 



1 



1 



1 



6irqaS(q) S{q) (2tt) 3 
1 - (k ■ qf 



d d k 



qk 2 



S(k + q)-l), 



(48) 



where we again have expressed the monomer diffusion 
constant D in terms of a Stokes radius a. After trans- 
forming to the dimensionless integration variable 



q 



q 

and performing the angular integration, one ha; 
1 



(49) 



C{q) 



6nqaS(q) 
1 1 

+ W)W) 2 Jo 



(50) 



dxf(x) (S(qx) - 1) 



with 



fix) = x 2 



1 + x 2 
2x 



hi 



1 +x 



l-x 



- 1 



This function can be expanded as 



/(*) = E 



i 



i 



2n + 1 2n + 3 



for x < 1, and 



/(*) = E 



n=0 



1 



1 



2n + 1 2n + 3 



(51) 



(52) 



(53) 



for x > 1. 

In order to make further progress, we have to spec- 
ify the structure factor S(q). This shall be done by the 
most simplistic model which takes into account both fi- 
nite bead size and finite chain length (see also Fig. [j]): 



N 



q < ^-N 



2tt at—v 



%N-»<q<^ 
9>¥ 



271 



(54) 



We now calculate C{q) in the scaling regime Rq 1 <C q <C 
a -1 . Defining the x values where S{qx) changes its be- 
havior as 



2tt 

Xl = frtF « l! 

2tt 

x 2 = — > 1, 
qa 



(47) we can write (qa)- 1 = x 2 /(2%), S^q)' 1 



N/S(q) = x~ 1/u ; hence 



(55) 
(56) 

-l/u 
c 2 i 



1 



, / dxf(x) 



W 1 JO 



(2^) 2 7 X1 



dxf(x)x~ 1/v 



1 f^ 2 
~ ttl^^ 1 " / dxf(x). 

( 27T r Jo 

Since x\ <C 1 and X2 3> 1, we can write 

cxi 4 

15 

da;/ (a;) » -x 2 , 

o 



(57) 

(58) 
(59) 



where we have taken just the leading-order terms of the 
expansions of /; this results in 



1 



15^ 2 Xl 



5-1/1/ 



(2^) 2 J X1 

1 1-1/1/ 
3^ 2 ' 



dxf{x) x - l/u 



(60) 



In the asymptotic limit (/-Rg — ► 00, i. e. 2,1 — ► 0, and 
qa — > 0, i. e. 2:2 — > 00, this obviously converges to the 
asymptotic value 



Gas — 



1 



dxf{x)x- 1 /" 



(27T) 2 Jo 

f 1/16 = 0.0625 i/=l/2 

~ \ v/3/(7tt) « 0.0788 v = 3/5 ' 1 j 

Focusing now on the correction, i. e. AC(q) = C(q) — C as , 
we find 



AC(g) 



1 



12tt 
1 



2^2 



l-l/i/ 



5-1/w 



+ 15^ ^ 

X /" x l 



(27T) 2 y 

1 



(27T) 2 y X2 

1 1-1/1- 

37r2 X 2 



dxf{x)x- 1/u 
dxf{x)x- 1/u 



(62) 
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taking again the leading-order terms for the remaining 
integrals results in 



AC(q) « - 



1 3 + v x-x/v 



12tt 2 1 
1 



1 



5-1/1/ 



15tt 2 5i/- 1 



(63) 



One thus sees that both finite chain length and finite bead 
size have the tendency to decrease C, as observed in the 
experiments. The latter effect is clearly more important, 
as xx enters only via a relatively high power. Further 
insight is gained by numerical evaluation of the shift for 
reasonable parameter. values. 

Tsunashima et alx3 performed their experiments with 
polyisoprene chains of size Rq = 210nm. Typical scat- 
tering wavenumbers in their plateau regime were given 
by qRc = 4 ... 8; the experimental observation in this 
regime was C w 0.06, i. e. a shift of AC «-2x 10~ 2 . 
In what follows, we consider the value qRa ~ 6. Thus 



x 2 = — w 500, 

qR G a 



(64) 



where xe have estimated the monomer size a as 
0.45nmt3. Inserting this into Eq. we find for the 
X2 contribution a value of AC « — 1 x 10~ 3 , i. e. one 
order of magnitude smaller than the experimental value. 

It is not completely clear if a more thorough treat- 
ment of the integral would fully account for the devia- 
tion; our guesSpia that it would probably not. Molecular 
Dynamics dataE3 seem to rather indicate that for typical 
systems (i. e. on not yet asymptotic length scales) the 
coupling between polymer and solvent is more complex 
than the simple Kirkwood description. Nevertheless, we 
consider our result as important for the interpretation 
of experimental data: There is obviously a substantial 
contribution to AC which stems from the finite bead 
size, and which is only weakly g-dependent. A plateau- 
like shape of C(g) alone apparently does not guarantee 
asymptotic behavior. Clearly more work has to be done 
to fully resolve the puzzle, but we believe our consider- 
ations show that theories which neglect the influence of 
finite bead size (and, to a lesser degree, of finite chain 
length) are simply not accurate enough to describe ex- 
perimental data even of rather long chains. 



III. NUMERICAL RESULTS 

In our numerical studies, we have used three different 
polymer chain models, which we will denote as model A, 
B, and C. 

Model A is a bead-spring model in the continuum. N 
monomers are connected via an anharmonic ( "finitely ex- 
tensible nonlinear elastic") spring potential, 



Ufene 



\kRl\xx 



Ro 



r<R 
r>Ro 



(65) 

where we use the standard parameters^ fc = 30, Ro = 1.5 
in dimensionless units. Between all monomers there is an 
additional non-bonded potential 

12 / -, \ 6 



-A, r < 2 1 / 6 



J COS 



|A [cos(ar 2 + /3) - l] 2 1 / 6 < r < 1.5 
0, r > 1.5, 



(66) 

where a and f3 are determined as the solutions of the 
linear set of equations 

2 1/3 a + /? = 7r (67) 
2.25a + /3 = 2tt, (68) 

i. e. a = 3.1730728678 and (3 = -0.85622864544. This 
potential has originally been constructed to simulate am- 
phiphilic systemsla. The parameter A serves to control 
the strength of the attractive interaction and is varied 
instead of the temperature, which is fixed at UbT = 1. 
For sufficiently strong A, the chain assumes a collapsed 
state, while A = corresponds to good .solvent. We used 
a combination of stochastic dynamicsEj and the pivot 
algorithm!! Applying standard methodscJ on data of 
chains of length of up to N = 2000, we located the 
point at A = 0.65 ± 0.02. In the good solvent limit, and 
at A = 0.65, we also ran an N = 5000 chain. 

Model B is a mesoscopic model for an aqueous solu- 
tion of the sodium salt of poly (acrylic acid) (PAA), 
whose input parameters have been derived from an ex- 
tensive atomistic simulation of an aqueous PAA solu- 
tion (T = 333.15 K and p = 1 atm) in the highly di- 
luted regime, such that the ion concentration (numher 
of charges on the chain, plus counterions) is 0.4 mol/lcil. 
From this simulation, structural averages like the distri- 
butions of bond angles or radial distribution functions 
between monomers were extracted. We mapped this sys- 
tem to the mesoscale by replacing one repeating unit (i. e. 
one monomer) by one bead. As center of the coarse- 
grained (CG) beads, the monomer center of mass (ex- 
cluding the sodium ion) was chosen. Bonded as well 
as non-bonded terms were parameterized by systemat- 
ically varying the interactions uptil the structure of the 
atomistic model was reproduced!! 2 ]. This also allowed us 
to neglect all explicit water molecules and sodium ions 
(necessarily present in the parent atomistic simulation) 
in subsequent CG simulations. Their effect on the PAA 
chain conformation is, however, implicitly present in the 
model. This means that a system of roughly 10 4 atoms 
could be reduced to a system which consists of only 23 
"super atoms" . As in model A, we used both stochastic 



8 




5000 



FIG. 2. R G = {R 2 g ) 1/2 and Rh = (i?^ 1 ) 1 for model A 
at G condition A = 0.65. Error bars are smaller than symbol 
sizes. Also shown are the fits as discussed in the text. 




0.15 



FIG. 4. Model A: Rg/Rh as a function of the scaling vari- 
able TV -1 / 2 , at O condition A = 0.65. The line results from 
the combined fits of Rg and Rh- 
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N 

FIG. 3. Model A: Same as Fig. [| but for good solvent 
condition A = 0. 

dynamics and pivot Monte Carlo moves. The final 
force field was utilized to calculate Rq, Rh and other 
static properties like the structure factor for PAA strands 
of length 8 to 3155 repeating unitaU. The numerical re- 
sults agree well with light scattering data on dilute PAA 
solutions with corresponding mean molar weights. In 
particular, the hydrodynamic radii of six different PAA- 
salt samples with molecular weights in the range from 
18100 to 296600 g/mol were measured. For four sam- 
ples, the molar masses M\y and the radii of gyration Rq 
were measured as well. The PAA samples were of poly- 
dispersity Dp between 1.5 and 1.8 and diktted in aqueous 
NaCl-containing solution (0.1 — 1 mol/l)0. 

Finally, model C is the SAW on the face-centered cubic 
lattice, which we prefer over the simple cubic for reasons 
of increased local flexibility, which in turn means prox- 
imity to the SAW fixed point. Units of length are defined 
in such a way that the bond length is v2- The chains of 
length = 64, 128, . . ..32768 were generated by using a 
dimerization procedureB. Up to N = 8192 the statistical 
sample always consisted of M — 1024 chains, while 




N -d-v, 

FIG. 5. Model A: Rg / Rh as a function of the scaling vari- 
able iV - ^"', at good solvent condition A = 0. The line 
results from the combined fits of Rg and Rh- 

M = 1085 for N = 16284 and M = 296 for N = 32768. 

In what follows, we outline our Rq and Rh data for 
these three models. Figures ^ and ^| summarize our re- 
sults for model A at 8 condition A = 0.65 (Fig. |^) 
and at good solvent condition A = (Fig. ||), respec- 
tively. For the chains, we obtained very good fits 
with the functions (R%) = 0.2834N _ 0.53 and (i?^ 1 ) 
2.71 Y 1 '' - 3.74N- 1 , while for the good solvent data 
the analogous fits are (R%) = 0.2706N 11754 - 0.32A^ 62 
and (R^) = 3.131N- - 5877 - 3.04N- 1 . These fit func- 
tions are also shown in Figs. § and |[ The ratio 
p = Rq/ Rh, as it results from these data, is shown in the 
subsequent Figs. || and || for and good solvent condi- 
tions, respectively. It should be noted that the numerical 
resolution (for each of our models) iSpClearly by far not 
competitive with the study by Li et alia. For this reason, 
we did not attempt to determine the exponents from our 
data, but rather used the values for v and A from Ref. 
|l|. We did not include an ]V _A term in the fit for Rh 
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FIG. 6. Model B: Scaling behavior of poly (acrylic acid) 
as measured by light scattering experiments and computer 
simulations with a coarse-grained model: (a) Radius of gy- 
ration Rg, (b) hydrodynamic radius Rh, (c) dimensionless 
ratio Rg/Rh- 

in the SAW case, although such a term is expected to 
be present. The reason is that our model A data are too 
inaccurate to allow for such a three-parameter fit in a 
stable way. Similarly, we ignored the non-analytic cor- 
rections to scaling in the O case, for essentially the same 
reason, as has been discussed in some more detail in Sec. 
[| Taking the statistical inaccuracies of the data, and of 
the resulting fit parameters into account, we obtain for 
the asymptotic amplitude ratio p = Rg/Rh the values 
p = 1.44 ±0.01 at the 6 point, and p = 1.63 ±0.01 in the 
excluded-volume case. The actual error in p is expected 
to be significantly larger, since neither the uncertainties 
in the exponents and in the location of the point, nor 
systematic errors due to higher-order corrections to scal- 
ing have been taken into account. This is particularly 
apparent in the O case, where one expects in the asymp- 
totic long-chain limit rather the Gaussian value 1.5045, 
but also obvious in the SAW case, where the results on 
the longer chains of model C yield a considerably smaller 
value for p. 

The most interesting aspect of model B is that it closely 
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FIG. 7. Model C: (R 2 G ) /{AN 2v ) as a function of iV" A , 
where we use A = 0.56, and A = 0.3341 from the (also shown) 
fit {R 2 G ) =AN x - mA +BN 0S2 . 

resembles a real system, and. a quantitative comparison 
with experiments is possibleO. In Fig. |^ we show simu- 
lation results for Rq, Rh, and their ratio. The data are 
taken as published in Ref. O. For the ratio, experimental 
results are also included. The scaling N^^ u \ and the 
extrapolation to p = 1.61 ± 0.02 is nicely borne out by 
the simulation data. The experiments are too inaccurate 
to demonstrate a clear systematic trend. In spite of this, 
an extrapolation yields p ~ 1.5 — 1.6, which means that 
the theoretical calculations are supported by data of a 
real chemical system. 

Our model C data (SAW) comprise the largest range 
of chain lengths of our three models, combined with pre- 
cise estimates of statistical errors, which allows a more 
detailed data analysis. For our Rg data, we obtained the 
fit (R 2 G ) = A7V 1 1754 + SiV a62 with A = 0.3341 ±0.0023, 
B = — 0.20±0.05, where we again use the exponents from 
Ref. [l]. The deviation \ 2 (sum of the residuals squares, 
normalized by the variances) has the value % 2 = 9.4 (10 
data points). The corresponding quality of fit Q, which 
is the probability to observe the measured % 2 value, or 
a larger one, is Q = 0.31. Our data, in a representation 
which emphasizes the corrections to scaling, are shown in 
Fig. [?]. It is seen that these are indeed weak, highlighting 
the difficulties in determining an accurate value for the 
correction-to-scaling exponent. 

Turning to our Rh data from model C, we first did a 
nonlinear two-parameter fit (Rjj 1 ) = AN~ v ' ft ' , result- 
ing in v e ff — 0.55. However, this fit is very poor, with 
a least-square sum \ 2 = 433. Conversely, a linear two- 
parameter fit (R^ 1 ) = AN^-^+BN- 1 yields a rather 
good value \ 2 = H-8 (Q = 0.16), with A = 2.732±0.005, 
B = — 3.10 ± 0.06, demonstrating also numerically that 
Rh data should be interpreted in terms of corrections to 
scaling, instead of an effective exponent. Actually, one 
should expect the presence of an additional correction of 
order iV~ A , A w 0.56. Since this correction tends to de- 
crease Rg (see Fig. [?]), it should also decrease Rh, i- e. 
increase (-R^} 1 ), or weaken the analytic N~( x ~ u } 



10 



O 0.4 




0.85 0.90 0.95 1.00 1.05 



FIG. 8. Model C: Quality of the fit (see text) for the rela- 
tion (Rh 1 ) = AN~°- 5Srr + BN~' t ' with fixed <j>, as a function 
of <b. 




0.20 



FIG. 9. Model C: (R^ 1 ) /{AN~ v ) as a function of 
iV"' 1 ""', where we use A = 2.753 from the fit 
(R^ 1 ) = AiV- ' 5877 + BN- 1 + CN- 1 - 15 , which is shown as 
well. 



term. Thus, in a regression (-R^ 1 ) = 0.5877 _)_ 

BN~^, where we keep <f> fixed, one should obtain the 
best fit for a value of 4> slightly smaller than unity. 
This is indeed what we observe, as seen from Fig. @, 
where we plot the quality Q of such a fit as a func- 
tion of <f). This figure also clearly rules out a sin- 
gle correction to scaling with an exponent of 1/2 or 
even larger. We thus attempted a three-parameter fit 
(i?^ 1 ) = AN- - 5877 ± BN- 1 + CN- 1 - 15 to also take the 
N~ A term into account. The result of this fit, which 
seems to be reasonably stable, is A = 2.753 ± 0.008, 
B = -4.3 ± 0.4, C = 2.2 ± 0.7, with X 2 = 5.0, and 
a very good quality Q — 0.66. We thus use this fit to 
demonstrate the corrections to scaling of {Rjj 1 } in Fig. 
|9|, where the presence of the N~ A term shows up in a 
slight curvature. Finally, we also used this fit, combined 
with the corresponding one for Rg (see Fig. |?]), to de- 
scribe the data on the ratio p = Rq/Rh, as shown in 
Fig. 10, where the asymptotic value is 1.591 ± 0.007. 
Again we feel that the real uncertainty is larger, due to 
lack of control of the systematic errors. We also checked 
that both the quality of fit, and the value of p did not 
change significantly pwhen we reduced the exponent A to 
its theoretical valueB A = 0.482. 

To summarize, we have collected our most important 
numerical results, the extrapolated p values, in Table ffl. 



Model 


Rg/Rh 


A (SAW) 


1.63 ±0.01 


B (SAW) 


1.61 ±0.02 


C (SAW) 


1.591 ±0.007 


A (6) 


1.44 ±0.01 



TABLE I. Asymptotic universal ratio Rg / Rh as estimated 
by numerical simulations of various models (see text). Error 
bars take into account statistical uncertainties only, while sys- 
tematic errors in the extrapolation procedure are neglected. 




FIG. 10. Model C: -Rg / Rh as a function of the scaling vari- 
able iV - ' 1- "' 1 . The line results from the combined fits of Rg 
and Rh- 
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APPENDIX A: EULER-MACLAURIN FORMULA 

Quite usually, sums are approximated via the <Wp% 
sponding integrals. The Euler-Maclaurin formulaElca, 
which we outline here for the convenience of the reader, 
constructs a systematic asymptotic expansion around 
that approximation. Defining a difference operator A 
via 



Af(x) = f(x + l)-f(x), 
one obviously has 

AF(N) = f(N) 

for 

JV-l 



F(N) = /(»)> 



and thus 



F(N) = A-\f(N) + const.. 
On the other hand, 

d 



A = exp 



dx 



1 



(Al) 



(A2) 



(A3) 



(A4) 



(A5) 



or 



A" 1 = 



dx 



exp 



dx 



-i -l 



- 1 



fe=0 



(A6) 



where Bk are the Bernoulli numbers defined via the Tay- 



lor expansion of x/(e x — 1): B a = 1, Bi 
B 2 = 1/6, B A = -1/30, . . . , B 3 = B 5 = B 7 = 
Hence, 



A" 1 



and thus 



dx 



1 1 d 



1 



12 dx 720 V dx 



-1/2, 
= 0. 



(A7) 



JV-l „N i 

f(n) = / dxf(x) - -f(N) + const. 

Jnn * 



1 d f{ \ 

12 Tx f{x) 



x=N 



720 dx* n ' 



x—N 

(A8) 



where the "integration" constant is determined via (per- 
haps numerical) comparison of both sides. For a power 
law with q < — 1 one thus finds from the definition of the 
Ricmann zcta function 



N-l 



N q+1 1 1 

Ni + a-q) + —qN 11 - 1 

q + 1 2 SV q > I2 q 



1 



^ q (q-l)(q-2)N^ 



(A9) 



By analytic contiimation with respect to q, this result 
holds for general cJ13. 



1 B. Li, N. Madras, and A. D. Sokal, J. Stat. Phys. 80, 661 
(1995). 

2 M. Doi and S. F. Edwards, The Theory of Polymer Dy- 
namics (Clarendon Press, Oxford, 1986). 

3 A. D. Sokal, in Monte Carlo and Molecular Dynamics Sim- 
ulations in Polymer Science, edited by K. Binder (Oxford 
University Press, New York, 1995), p. 47. 

4 G. Besold, H. Guo, and M. J. Zuckermann, J. Polym. Sci., 
Polym. Phys. 38, 1053 (2000). 

5 L. Schafer, Excluded Volume Effects in Polymer Solu- 
tions as Explained by the Renormalization Group (Springer, 
Berlin, 1999). 

6 J. Zinn- Justin, Quantum Field Theory and Critical Phe- 
nomena, 2nd ed. (Clarendon Press, Oxford, 1993). 

7 M. Adam and M. Delsanti, J. Physique 37, 1045 (1976). 

8 M. Adam and M. Delsanti, Macromolecules 10, 1229 
(1977). 

9 N. Nemoto, Y. Makita, Y. Tsunashima, and M. Kurata, 
Macromolecules 17, 425 (1984). 

10 K. Venkataswamy, A. M. Jamieson, and R. G. Petschek, 
Macromolecules 19, 124 (1986). 

11 D. Reith, B. Miiller, F. Muller-Plathe, and S. Wie- 
gand, preprint, submitted to Macromolecules, cond 



mat/0103329 . 

Y. Tsunashima, M. Hirata, N. Nemoto, and M. Kurata, 
Macromolecules 20, 1992 (1987). 

13 R. M. Jendrejack, M. D. Graham, and J. J. de Pablo, J. 
Chem. Phys. 113, 2894 (2000). 

14 G. Weill and J. des Cloizeaux, J. Physique 40, 99 (1979). 

15 L. Schafer and A. Baumgartner, J. Physique 47, 1431 
(1986). 

16 K. Kremer, A. Baumgartner, and K. Binder, Z. Phys. B 
40, 33 (1981). 

17 Y. Oono and M. Kohmoto, J. Chem. Phys. 78, 520 (1983). 

18 Y. Tsunashima, Macromolecules 21, 2575 (1988). 



12 



19 J. F. Douglas and K. F. Freed, Macromolecules 17, 2354 

(1984) . 

20 C. M. Guttman, F. L. McCrackin, and C. C. Han, Macro- 
molecules 15, 1205 (1982). 

21 C. M. Guttman, J. Stat. Phys. 36, 717 (1984). 

22 A. Z. Akcasu and C. M. Guttman, Macromolecules 18, 938 

(1985) . 

23 J. Hager and L. Schafer, Phys. Rev. E 60, 2071 (1999). 

24 P. Grassberger and R. Hegger, J. Chem. Phys. 102, 6881 
(1995). 

25 M. Schmidt and W. Burchard, Macromolecules 14, 211 
(1981). 

26 Y. Tsunashima et al, Macromolecules 20, 2862 (1987). 

27 J. Batoulis and K. Kremer, Macromolecules 22, 4277 
(1989). 

28 A. J. C. Ladd and D. Frenkel, Macromolecules 25, 3435 
(1992). 

29 A. Z. Akcasu and H. Gurol, J. Polym. Sci., Polym. Phys. 
14, 1 (1976). 

30 M. Benmouna and A. Z. Akcasu, Macromolecules 13, 409 

(1980) . 

31 C. C. Han and A. Z. Akcasu, Macromolecules 14, 1080 

(1981) . 

32 M. Bhatt and A. M. Jamieson, Macromolecules 21, 3015 



(1988). 

33 M. Bhatt, A. M. Jamieson, and R. G. Petschek, Macro- 
molecules 22, 1374 (1989). 

34 Y. Shiwa, J. Phys. A, Math. Gen. 24, L579 (1991). 

35 Y. Tsunashima, Polym. Journ. 24, 433 (1992). 

36 B. Diinweg and K. Kremer, J. Chem. Phys. 99, 6983 
(1993). 

37 B. Hammouda and A. Z. Akcasu, Macromolecules 16, 1852 



(1983). 

Y. Tsunashima, private communication. 



39 K. Kremer, in Monte Carlo and Molecular Dynamics of 
Condensed Matter Systems, edited by K. Binder and G. 
Ciccotti (Italian Physical Society, Bologna, 1996), p. 669. 

40 T. Soddemann, B. Diinweg, and K. Kremer, Europ. Phys. 
J. E 6, 409 (2001). 

41 O. Biermann et al. } preprint, submitted to Macromolecules, 



cond-mat/0101115 



42 D. Reith, H. Meyer, and F. Miiller-Plathe, Macromolecules 
34, 2335 (2001). 

43 F. W. J. Olver, Asymptotics and Special Functions (Aca- 
demic Press, New York, 1974). 

44 R. L. Graham, D. E. Knuth, and O. Patashnik, Concrete 
Mathematics (Addison- Wesley, Boston, 1994). 



13 



